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We present a simple analytical method to solve master equations for finite temperatures and any 
initial conditions, which consists in the expansion of the density operator into normal modes. These 
modes and the expansion coefficients are obtained algebraically by using ladder superoperators. 
This algebraic technique is successful in cases in which the Liouville superoperator is quadratic in 
the creation and annihilation operators. 
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In general, any physical system is affected by its sur- 
roundings since, even in laboratory, it is hard to avoid 
its interaction with the environment. Furthermore the 
dissipative processes and the continuous fluctuations of 
the environment generate decoherence processes, which 
are responsible to an irreversible loss of information of 
the system. Thus the environment plays a fundamental 
role in the dynamics of any system so that in any real 
physical situation it is necessary to analyze its effects. 

In a Hamiltonian formalism, certain types of reservoir 
may be simulated by a bath of harmonic oscillators lin- 
early coupled to the system and, in some cases, the solu- 
tion of the full problem of coupled system-reservoir dy- 
namics may be exactly obtained in the density operator 
formalism [1]. Nevertheless, we are in general only inter- 
ested in observing the system so that we should eliminate 
the environment by tracing over the bath variables, which 
provides us the reduced density operator of the system, 
p{t). 

The approximated dynamics of the reduced density op- 
erator does not need of the full solution: we may gener- 
ate, through the trace over the bath and suitable approx- 
imations, master equations that describe its evolution. 
These equations become irreversible due to the fact that 
the system-reservoir correlations are suppressed in the 
reduced density operator. The most common approx- 
imation consists in adopting the Born-Markov limit in 
which the system-bath coupling is weak and the mem- 
ory of the reservoir is short. In this limit, the master 
equation becomes local and the effects of the reservoir 
are resumed in the form of a non-Hamiltonian Lindblad 
term. The general form of this type of equation is [2,3] 



dp{t) 
dt 
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(1) 



where H is the Hamiltonian of the system and Cp is the 
Lindblad term generated in the elimination of the bath 
degrees of freedom. 

Several analytical methods to solve master equations 
of this type for any initial condition can be found in 
Refs. [4-10]. The simple case of an oscillator in a 
zero-temperature reservoir has been solved by Guerra et 
al. [4] . Mokarzel [5] has solved the case of two coupled 
oscillators at zero temperature and in the presence of 
pump in both modes by using a method based on the 
Lie algebra for superoperators. This method consists in 



decomposing the exponential of the total Liouvillian into 
simpler exponentials. Gilles and Knight [6] have solved 
a master equation in which the Lindblad form contains 
terms non-quadratic in the creation and annihilation op- 
erators, like a?pa^ ^. A wide class of problems at zero 
temperature, including Hamiltonians and Lindblad terms 
non-quadratic in the operators a and a\ has been solved 
by Klimov and Romero [7]. In this case, superopera- 
tors belonging to a deformed algebra have been used. 
Walls and Milburn [8] have used quasidistributions to 
obtain the solution of an oscillator coupled to a thermal 
reservoir without the usual rotating- wave approximation. 
Pefinova and Luks [9] also have used quasidistributions 
to solve a master equation at finite temperature with a 
Hamiltonian of the form ^a^. Briegel and Englert [10] 
have elaborated a method based on the expansion of the 
density operator of the problem into the eigenstates of the 
Liouvillian. The eigenstates have been obtained by us- 
ing differential equations and, the expansion coefficients, 
by taking a scalar product of the initial density opera- 
tor with a dual basis orthonormal to the eigenstate basis. 
This method has been used [10] to solve some problems 
as a single oscillator in a finite temperature bath and a 
two-level system coupled to a cavity field mode through a 
Jaynes-Cummings type interaction at zero temperature. 

In this paper we present an algebraic technique to ob- 
tain the solution of master equations, which consists in 
the expansion of the density operator into the eigenstates 
of the Liouvillian. The eigenstates and the expansion co- 
efficients are easily obtained in an algebraic way: we de- 
fine step-down superoperators, which are used for obtain- 
ing the steady state and the expansion coefficients, and 
step-up superoperators, which enable us to obtain the re- 
maining eigenstates. In the next sections we present some 
examples where this technique is successful: harmonic os- 
cillator in a zero-temperature bath and in a thermal bath, 
two-level system in a thermal bath and two coupled os- 
cillators in a zero-temperature bath and in a finite tem- 
perature bath. To our knowledge the last case has not 
yet been solved. 



I. DAMPED HARMONIC OSCILLATOR 



Eq. (1) may be rewritten in a more compact form 
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Using the above relations and the definition of /C, we 
may obtain the commutation relations 



where the Liouvillian superoperator K is defined as 

1C=^{H.- .H)+C . (3) 
tn 

Eq. (2) has the formal solution, 

p{t)=e'^*p{0). (4) 

However, it is not trivial, in general, to obtain a useful 
expression for the density operator. To allow the action 
of this time evolution operator wc will expand the initial 
density operator into the eigenstates of the Liouvillian. 

In order to exemplify our method, we will analyze the 
behavior of a harmonic oscillator in contact with a zero- 
temperature reservoir. Then, we will solve a master equa- 
tion of type (1), where the Hamiltonian of the system is 
given by 

H = hioa^a (5) 

and the Lindblad superoperator that simulates the reser- 
voir has the form [2,3] 

jC = ^{2a..a^ - a. - .a^ a) , (6) 

being 7 the damping rate of the system energy due to 
the action of the reservoir. 

Here we have adopted the notation A. {.A) for super- 
operators that represent the simple action of an opera- 
tor A to the left (right) on the target operator, p, i.e., 
A.p := Ap {.Ap := pA). In order to avoid confusion be- 
tween operators and superoperators, the latter will be 
denoted by either a Calligraphic letter or a Roman letter 
followed or preceded by a point. It is easy to demonstrate 
that superoperators of this form obey the following rela- 
tions: 

A.B. = AB. , .A.B = .BA, 

.BA. = A..B , .1 = 1. = 1 , (7) 

which, obviously, only have a meaning if acted on an 
arbitrary target operator. 

By defining the commutator between two superopera- 
tors in the usual form, [A^ B] :— AB — BA, we can show 
that almost all the commutator properties among oper- 
ators remain valid here, for example, 

[\A + 5B,C] = \[A,C] + 5[B,C], 

[B,A] = -[A,B], 
[AB,C]^[A,C]B + A[B,C]: (8) 

where A and 6 are complex numbers. Furthermore, from 
relations (7), it is easy to obtain the following properties: 

[A.,B] = [A,B]. , [.A,.B] = .[B,A], 
[A.,.B]=0 . (9) 



[/C, aK] = {-iuj - 7/2) + 7 .a'f , 
[lC,.a^] = (-iw + 7/2) .a' , 

[/C, .a] = [iiv — 7/2) .a + j a. , 

[]C,a.] = {iu + -f/2) a. . (10) 

We may obtain from linear combinations of these rela- 
tions, new superoperators, 

M+=a^.-.a^ , M-=a., 

Af+ = .a-a. , Af-= .a^ (11) 

which satisfy simpler commutation relations with /C, 

[/C,A^±] = ±(-ia;-7/2)A^±, 
[}C,M±] = ±{iuj - 7/2);V± . (12) 

Therefore M.± and A/± are ladder superoperators. The 
definition of these ladder superoperators was possible be- 
cause the commutation relations (10) correspond to a 
closed algebra. In fact, this technique is limited to cases 
in which the commutation relations between the Liou- 
villian and the superoperators that constitute it form a 
closed algebra. 

The ladder superoperators were normalized so that 
they satisfy commutation relations identical to those of 
two decoupled harmonic oscillators, 

[M.,M+] = 1 , [A/-_,A/-+] = 1, 

[M-,M+] = , [M-,M-]=0. (13) 

Using these superoperators we may rewrite the Liou- 
villian JC in the suggestive form, 

/C = {-iuj - 'y/2)M+M- + {iuj - j/2)J\f+J\f_ . (14) 

Due to the fact that and M+N- commute, 

we can find simultaneous eigenstates of these superoper- 
ators, i.e., 

M+M- i?™'" = m i?"'" , 
J\f+J\f- i?™'" = n i?™'" . (15) 

It is clear that these eigenstates will be also eigenstates 
of the Liouvillian, 

K BJ^''' = [m{-iijO - 7/2) + n{iuj - 7/2)]i?"'" . (16) 

Since the density operator determines the ensemble 
average of observables by means of traces, {A{t)) = 
tr(Ap(f)), we must assume that it has finite trace with 
any other operator. More generally, the density operator 
should have finite trace with any superoperator. Since 
the eigenstates i?™'" will be used in the expansion of the 
density operator, we must impose that they also have this 
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property. Therefore, solutions of Eq. (15) like 1 (with 
eigenvalues m = — 1 and n — —1), a [m = —2,n = —1), 
(m = — 1, n = —2) and aa^ (to = —2, n = —2) cannot 
represent i?'"'"'s, because they have divergent trace with, 
for example, the superoperators 1., a^., a. and 1., respec- 
tively. This shows that negative eigenvalues do generate 
unphysical solutions. In fact, these solutions correspond 
to density operators that grow without bound as t in- 
creases, which is not physically expected in the case of 
a zero-temperature reservoir. Therefore, we will only at- 
tempt solutions of Eq. (15) with positive eigenvalues. 

From the commutation relations (13) and from 
Eq. (15) we can show that A1+ increases the eigenvalue 
of Mjf^M- by one unity and does not change the eigen- 
value ofAf+Af-. In analogous form we have that A/V only 
changes the eigenvalue of Af+Af- . So far we have not spec- 
ify any type of normalization for the eigenstates i?*"'". 
Therefore, we may choose the coefficients of the step-up 
relations as being \/m + 1 and ^/rT+T, in analogy with 
the harmonic oscillator, so that the step- up relations have 
the form 



min(m,n) 



Af+ ff"'" = v^m . 



(17) 



From relations (13), (15) and (17) we obtain the step- 
down relations 



AT- iZ'"'" = i?"'"-! . 



(18) 



We will look for the eigenstate of null eigenvalues, i.e., 
R°'^. The above relations show that we cannot obtain 
new eigenstates by applying A4- and jV- to BP'^, 

M- = aijO'O = , ^_ = ijO'Oat = . (19) 

Although the convention used in the step-up relations 

gives a restriction in the normalization of the eigenstates, 
we still have the freedom to choose the normalization of 
one of them, say R^'^. We will choose it such that 

tri?°'0 = 1 . (20) 

Thus, this eigenstate takes the simple form 

i?°'° = |0)(0|, (21) 

where |0) represents the oscillator ground state. The 
remaining eigenstates with positive integer eigenvalues, 
m,n = 0, 1, 2, ... , may be then obtained by applying 
A4+ and Af+ successively to 



! \/m\ 



(22) 



Thus, we may obtain the explicit form of the eigenstates 
in the Fock basis 



fe=0 



(-1)'= 



7 I 



(to — fe)! (n — fc)! 
x\m-k){n-k\. (23) 



Using the cyclic property of the trace, we get 
tr(A^+ p) = and tr(A+ p) = for any target opera- 
tor p. Thus, taking the trace of Eq. (15), we obtain that 
TO tri?™'" = and n tri?"*'" = 0. This allows us to con- 
clude that the trace of i?™'" will be identically null for 
TO or n different from zero. Taking into account the con- 
vention (20), we may write 

tri?™'" = 6m,0 Sn,o ■ (24) 

By decreasing the eigenvalues of a given eigenstate, 

Afi Ml ' 



R" 



ml 



i\{m - iy. j\{n - j) 



R" 



, (25) 



and taking the trace of this equation, we obtain a very 
useful property. 



tr 



Afi M'_ 



(26) 



It is shown in the appendix that the R"^''"'s with pos- 
itive integer eigenvalues, to,, n = 0, 1. 2, . . . , form a com- 
plete set, which allows us to expand any density operator 
into this eigenstates. In particular we may expand p(0). 



p(0) = X] • 



(27) 



m=0 n=0 



Using Eq. (26), we easily obtain the expansion coeffi- 
cients, 



/TV" M""- 
V vTi! vto! 



P(0) 



or more explicitly. 



tr(a'"p(0)a^ ") . 



(28) 



(29) 



The solution of the master equation will be obtained 
by substituting Eq. (27) into Eq. (4) and using Eq. (16), 



X e 



m— n— 

'm{—iio — 7/2) -|- n{iuj — 7/2)]t 



(30) 



Thus, to obtain the solution p{t) of a given initial con- 
dition, p{Q), we should determine the coefficients C""'" 
and perform the summations over the indices to and n in 
Eq. (30), using the explicit form of the eigenstates i?™'". 

Notice that C°'° = trp(O) = 1, independently of the 
initial condition. Thus, from expansion (30), it is clear 
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that any initial state decays to the eigenstate 7?°'" as t 
goes to infinity. This is due to the fact that the system 
dissipates energy to the reservoir until the system reaches 
its ground state. 

As an example, we will use as initial condition the os- 
cillator in a Fock state, 



P;v(0) = |Ar)(Ar|, 
where the expansion coefficients are given by 

AT! 



'"n!(7V-n)! ' 



(31) 



(32) 



By performing the summations in Eq. (30), we will obtain 
the density operator at time t, 



N 



k=0 



k\{N-k)\ 



T^)^-'=|fc)(fc|. 



(33) 



We can also easily calculate the density operator when 
the oscillator is initially in a coherent state, 

Pa(0) = |a)(a|. (34) 

In this case, we have for the expansion coefficients 



Vm! vn! 

and then the density operator will be given by 



\ae 



(-iw-7/2)t 



){ae 



(-iw-7/2)t| 



(35) 



(36) 



indicating that the oscillator dissipates but remains as a 
pure state during all the evolution. 



II. HARMONIC OSCILLATOR IN A THERMAL 
BATH 

We now consider the case of a harmonic oscillator in 
contact with a thermal reservoir. The Hamiltonian of the 
system is the same as that given by (5), but the Lindblad 
superoperator is now given by [2,3] 



£ = (n + l)^(2a..o^ - a^a. - .a^a) 
n^{2a^ ..a — aaK — .aa^) , 



(37) 



where n = (e ~ 1)^^ is the average number of quanta 
in the reservoir, /3 = l/fc^T, being ks the Boltzmann's 
constant. For this case, we will define new step-up and 
step-down superoperators, 

A4+ = a^. — .a^ , A4- = — ft .a + {ft + 1) a. , 

Kf+ = .a-a. , Kf- = -na^. + {n + 1) .a) . (38) 



Using these superoperators the new Liouvillian /C may 
be rewritten in a form analogous to Eq.(14), 

K. = {-iuj - '^/2)M+M_ + {iuj - '^I2)N+Kf- . (39) 

As the superoperators M+, A4-, A/+, and Af- obey 
the same commutation rules as those obeyed by the 
superoperators M+, A4-, A/+, and 7V_, we follow the 
same procedure used in the last section to obtain the 
simultaneous eigenstates of M+M- and Af+Xf-. We 
denote these eigenstates by i?'"'"(n). The eigenstate 
i?^'''(n) must have unit trace and obey the equations 
M- R°'°{n) = and A/L R°'°{n) = 0, which will give 
us 



fc=0 



(n + 1) 



fc+i 



(40) 



This operator exactly represents a thermal distribu- 
tion, i?0'0(n) = e-^-^/Z, being Z = tre"^^ the par- 
tition function of the problem, since this state is in 
thermal equilibrium with the reservoir. The remain- 
ing eigenstates with positive integer eigenvalues, m,n = 
0,1,2,..., may be obtained by applying A4+ and A/+ 
successively to i?°'°(n), 




Kfn kAm 



n! vmi 



(41) 



\kl 



k=0 
w p k,k—n 



m!(fc + TO-n)! (n -I- l)™+i 



n + 1 



\k + m — n){k\ , m> n 



(42) 



E 



m\kl 



1 



fc=0 

w p k,k—m 
n 



nl{k + n-my. (n -|- l)"+i 
n 



n + 1 



\k){k + n — m\ , m <n, 



where P ^ (x) are the polynomials 

pk,i(^^_ V (-iv-i—il±JI!:}—^ (43) 

j=max(0,l) 

Following the procedure adopted in the last section, 
the solution of the master equation for the case of a har- 
monic oscillator in a thermal bath takes the form 



/^T(i) = EE^™'"(^)^"'"(^) 

m=0 n=0 

,[m{-iu) - 7/2) + n{iLO - ^/2)]t 



X 



(44) 
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where the expansion coefficients are obtained in a form 
analogous to Eq.(28), 

C-."(n) =tr^^p(0) . (45) 

As a simple example, we consider an oscillator which 
is initially in a thermal state at temperature Tq different 
from the reservoir temperature T, 

MO)=E (^, + l).+i |fe)(fe|' (46) 

where no = (e^"^'^ - being /?o = l/ksTo. In this 

case, the expansion coefficients take the very simple form, 

C™'" = ^„,„(no-n)". (47) 

The solution of the master equation is obtained substi- 
tuting the above result into (44), rearranging the sum- 
mations and using the identity [11] 

So that the density operator remains as a thermal distri- 
bution during all the time evolution, 

m*) = E7^(|Vt)^i^)(^I' (49) 

where n{t) = n+{no—n)e~^^ , showing that the oscillator 
reaches the thermal equilibrium with the reservoir in a 
characteristic time I/7. 

III. TWO-LEVEL SYSTEM IN A THERMAL 
BATH 

Another interesting example is a two-level system, with 

a ground state \g) and a excited state |e), in contact with 
a thermal reservoir. In this case, we must solve a master 
equation of type (1), where the Hamiltonian of the system 
is given by 

H = ?iuja^/2 , (50) 

being cr^ = |e)(e| — |5)(fif|, andtheLindbladsuperoperator 
is given by [2,3] 

- T 

jC = {n+ 1) — (2(T_..(T+ — (T+(T_. — .(T+(T_) 

-|- n— (2cr_|-..cr_ — cr_cr_|_. — .(T_(T+) , (51) 

where cr+ = \e){g\, (7_ = \g){e\ and n = (ef^^ - 1)^^. 
The master equation for this problem can be easily solved 
in the basis of the states |e) and \g). Although we are 



aware of this, we proceed to solve this problem with our 
method to give a simple example of dealing with anti- 
commutation relations. 

The Lindblad superoperator can be rewritten in a more 
convenient form, 

_ r 

£ = (1 - N)-{2a-..a+ - a+a-. - .cr+(7_) 

- r 

+ N-{2a+..a--a-a+.- .a-a+), (52) 

where N = {e>^^^ + l)'^ and T = 7 coth{f3nu/2). For 
this case, we will define the step-up and step-down su- 
peroperators 

V+ = a+. + Uz..u+ , 

V- = {1- N)(j_. + Na^..(7_ , 

Q+ = -(J- + <J--(Jz , 

Q_ = {l-N).a++Na+..a, . (53) 

By defining the anticommutator between superopera- 
tors in the usual form, {A, B} := AB + BA, we will have 
the following anticommutation relations: 

{V-,V+} = 1 , {Q_,Q+} = 1, 
{V+,V+} = , {Q+,Q+} = 0, 
{P_,P_} = , {Q_,Q_} = 0, (54) 

and the following commutation relations: 

[P+,Q+] = , [P+,Q_] = 0, 

[V-,Q+] = , [P-,Q-]=0, (55) 

where 1 = |e)(e| -|- We may again rewrite /C in 

terms of these superoperators, 

/C = {-iuj - T/2)V+V- + [iuj - r/2)Q+Q_ , (56) 

and find simultaneous eigenstates of V+V- and Q+Q_, 
due to the fact that they commute, 

V+V- SP'"^ = p SP'" , 

Q+Q- SP'i = q SP'^ . (57) 

The anticommutation relations enable us to obtain the 
commutation relation [V+V-,V-] = —V-. Applying this 
relation on the eigenstate 5^'^ and using V^ = 0, we will 
have (p - 1)V-SP''^ = 0. Then, acting V+ on this last 
expression, we will obtain p{p — 1)3^''^ = 0, allowing us 
to conclude that, if we want 5^''' 7^ 0, we must have 
p = 0, 1. We may proceed in analogous form and obtain 
q = OA. 

Here wc also have step-up and step-down relations, 
which we will not explicit. The eigenstate S°'° is obtained 
from equations V-S°'° = 0, Q-S°'° = and tTS°'° = 1. 
Then we find 

S°'' = N\e){e\ + {l-N)\g){g\, (58) 
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which also can be identified as a thermal distribution 
S°'° = e-^^^/Z, where H = huj(j^/2 and Z = tre'^^ 
is the partition function of the problem. The remaining 
eigenstates will be obtained by stepping S^'^ up, 



which will give us 

S''^ = \e){g\ , S'-' = \9){e\, 
S'^' = \e){e\-\g){g\. 



(59) 



(60) 



Obviously, the basis formed by four eigenstates S^'' 
is complete since, in this case, the space of operators is 
four-dimensional . 

Following the procedure adopted in the first section, 
the solution of the new master equation takes the form 

p{t) =J2Y1 ^"'^ 'S^-'etPC-^^ - r/2) + - T/2)]t ^ 



p=0 g=0 



(61) 



where the expansion coefficients are obtained from ex- 
pression 



= tr(QlP^ p(0)) . 



(62) 



It is interesting to note that there is more than one 
choice for ladder superoperators. For example, the su- 
peroperators 

P+ = .cr+ - (T+..az , 
V'_=N.(T_ - (l-7V)o-_..CT^, 



Q'_ = Na+. - {1 - N)a,..a+ . 



(63) 



obey anticommutation relations identical to (54) and 
commutation relations identical to (55). Furthermore, 
these superoperators rewrite the Liouvillian in a form 
which is analogous to (56), due to equalities V'^V_ = 
V+V- and Q'^Q'_ = Q+Q-. However, the superopera- 
tors P+, V-, and Q'_ have anticommutation relations 
identical to (54) but, instead of obeying the commuta- 
tion relations given in (55), they obey anticommutation 
relations. An analogous fact occurs with the group of 
superoperators "P:^, V'_^ Q+ and Q_. 

As a simple example, we can use again as initial con- 
dition the system in a thermal state at temperature Tq 
different from the reservoir temperature. 



p{Q) = No\e){e\ + {l-No)\g){g\, 



(64) 



where Nq = {e^°^'^ -|- 1) i. In this case, we will have for 
expansion coefficients 

Cf'« = 5p,g(iVo-7V)«, (65) 

so that the density operator at time t will be given by 



p{t) = N{t)\e){e\ + {l-N{t))\g){g\, 
-Tt 



(66) 



where N{t) = N + {No - N)e i.e., the system re- 
mains as a thermal distribution during all the evolution, 
reaching the thermal equilibrium with the reservoir in 
the characteristic time l/F. 



IV. TWO COUPLED MODES IN CONTACT 
WITH A ZERO-TEMPERATURE RESERVOIR 

Wc will analyze now the case of two (;ouplcd harmonic 
oscillators in contact with a zero-temperature reservoir. 
In this case, the master equation is still of type (1), but 
now the Hamiltonian of the system is given by 

H = tiUaO^a + fiubb^b + hga^h + ng*ab^ (67) 

and the Lindblad of this problem is 

" a'a. — .a'a) 



C = -f{2a..a^ 



+ ^(26..6t-6t6.-.6t6) 
+ ^(2a..6l'-6^a.- .feto). 



(68) 



where \jc\ = \Jla)h for the case in which the reservoir 
is common to the two modes and 7c = for the case in 
which there are separate reservoirs for each mode [5] . For 
both cases we will define a set of step-up superoperators 

M+ = r+Ml + s+Ml , 
Ar^_ = r^M% + S-M^ , 
V+ = r\Ul + sXNX , 
Q+ = r*_Nl + s*_NX , 



(69) 



and step-down superoperators 

M- = u+Mt + v+M''_ , 
M- = u-M''_ + v^Mt , 
V- = u\Nl + v^Mt , 



(70) 



where the superoperators with indices a, b refer to those 
of a single mode of the form (11) (for example, M.°L = 
a. and M''_ = b.). The coefficients in (69) and (70) are 
defined as 



r± 



S±A 



u± = ± 



^2iA{S±A) 
S±A 



s± = i- 



V 



^/2iA{S±A) 



y'2iA{S±A) 

v± = 71 
^2iA{S±A) 
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where the parameters C/, S and A are given by 



cx;' oo oo oo 



U = g 

UJa - UJb 



,1a -lb 



2 4 



(72) 



These superoperators obey the commutation relations 
[M-,M+] = '^ , W-,M+] = l, 



[V-,V+] = 1 , [Q_,Q+] = 1, 



(73) 



whereas the remaining relations are identically null. The 
Liouvillian /C can be rewritten in terms of these super- 
operators as 



IC = X+M+M- + X-Af+Af- 
+ XIV+V- + X*_Q+Q- , 



(74) 



where 



. i^-^ . (75) 



Following the procedure used before, we may find the 
simultaneous eigenstates oi M+M-, Af+Af-, V+P- and 
Q+Q- with eigenvalues m, n, p and q, respectively, which 
will be denoted by _R™^"'P'9. Again, we have step- up and 
step-down equations, which we will not write explicitly, 
and ijO'O'O'O is the unit trace eigenstate whose indices 
cannot be decreased. Its explicit form is 



^0,0,0,0^^0,0^0,0^ 



(76) 



where = |0)o- o-(0|, with a = a,b, is the eigen- 

state of lowest eigenvalue for a single mode, according 
to Eq. (21). The remaining eigenstates are obtained by 
stepping RO'Ofifi 



nl -pP \fn KAm 
j^m,n,p,q _ ^+ ■'^ + ■'^'■+ j^O.0.0.0 

and their explicit forms will be 

m+n p+Q 



where 



(77) 



(78) 



k=0 1=0 



fc!(m + n — k)\ 
m\n\ 

n\ 



min(m,fe) 

E 

J— max(0,A;— n) 



j!(m- j)! 



{k - j)\{n - k + j)\ 



i m—j k—j n — k-\-j 

-r-!^s_^ •'r_ s_ 



(79) 



and i?™'", with a = a,b, are the eigenstates for the case 
of a single mode given in Eq. (23). Thus, the solution of 
the master equation for the case of two coupled modes in 
contact with a zero-temperature reservoir is 



m— n— p— q—0 

^ g(mA+ +nX- -|-pA+ + qX*_)t 



where the expansion coefBcients are 
/ Ql AfH 



: tr 



P(0) 



(80) 



(81) 



By defining the operators x± = u±a + v±b, we may 
rewrite these coefficients in a simpler form, 



^m,n,p,q 



1 



\/m\n\p\q\ 



tr(x!^x!r P(0) xtV-") • (82) 



As example, we may consider the case in which both 
modes are in a coherent state 

P(0) = \a)aa{a\\p)bb{f3\. (83) 
In this case, the expansion coefficients will be given by 



{uXa*+vlP*)P {u*_a* + v*_l3*)'i 



and the solution of the master equation is 

p{t) = \a{t))a a{a{t)\ \m)b b{P{t)\ , 

where 

a{t) = aF{t) + (3G{t) , 
I3{t) = aH{t) + pi{t) , 



(84) 



(85) 



(86) 



being 

(Fit) Git) 



S+ S- 





u+ v+ 



(87) 



With the coefficients (71), we will obtain the explicit form 
of these functions 



F{t) 
G{t) 



cos{At) - i^sen{At) ) e , 



-^— sen(At)e 



S_ 
A' 

-iRt 



H{t) = -i^sen{At)e-'^* . 



m 



S \ 
cos{At) + i— sen(Af) )e~*-^* 



(88) 



We see that, in this example, the two modes are always 
separable. 
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V. TWO COUPLED MODES IN CONTACT WITH resulting in 
A SINGLE THERMAL RESERVOIR 

Wc want now to include temperature in the system 
formed by two coupled modes in contact with a single 
thermal reservoir. It only will affect the Lindblad super- 
operator that will be now [5] 

C = {n + l)-^{2a..a^ — a^a. — .a^a) 
+ n^(2a^.o — aa^ . — .aa^) 
+ (n + l)y(26..6l'-6'f&. -.6'f6) 
+ n^(26^.6-66^-.66'f) 
+ {n + l)^(26..at - a^b. - .a^b) 
+ n^{2a^..b-ba^.- .ba'') 
+ (fi + l)^(2a..6t-6ta. -.6ta) 
+ n'^{2bl.a-ab''.- .ab''), (89) 



where |7c| = sj'jalb- The case of two reservoirs at the 

same temperature may be obtained taking 7c = 0. 
The step-up superopcrators are now given by 

M+ - r+Ml + s+M\. , 

V+ = rXNl + s\K[X , 

Q+ = r*_Kfl + s*_NX , (90) 

whereas the step-down superopcrators are 
M- = u+Mt + v+Mt , 



Q- 



(91) 



where the superoperators with indices a, b refer to those 
for the case of a single thermal mode of the form (38) 
(for example, Mt. — —n.a + (n + l)a. and M.^_ = —n.b + 
{n + 1)6.). The steady state for the case of two coupled 
thermal oscillators is factorable into the steady states of 
each mode of the form (40), 



i?0'0'0'0(n) = i?r(n)flr(n) 



(92) 



indicating that the steady state is the one in which each 
oscillator is in thermal equilibrium with the reservoir. 
The remaining eigenstates are obtained by stepping the 
steady state up, 

n? -pP A?" /Cf™ 
.9(n) = 4± 1± l^ijO,o,o,0(^) ^ (93) 

vQ- VP- V n! Vm! 



m+n p+q 

fe=o ;=o 

X El'J (n)ii™+"-'='^+«-' (n) , (94) 

where D™'" are given in Eq. (79) and -R^'"(n), with a = 
a, b, are the eigenstates of a single thermal mode given in 
Eq. (42). 

The expansion of the density operator into the eigen- 
states will be similar to Eq. (80), 

oo oo oo oo 

PT(t) = ^ ^^^C'"'"'P'«(n)i?'"'"'P'«(n) 

m— n— p—0 q—0 

^ g(mA+ + nX- + pX'^ + qX*_)t 
where the expansion coefficients are given by 



(95) 



(^) = tr(^-^^^p(0) 



n! vm! 



VI. CONCLUSIONS 



(96) 



To summarize, we have developed a simple algebraic 
method to obtain the solution of master equations with 
quadratic Liouville superoperators for finite tempera- 
tures and any initial conditions. This method consists 
in the expansion of the initial density operator into 
the eigenstates of the Liouvillian, allowing its evolution 
in a simple form. For obtaining these eigenstates we 
have used ladder superoperators, which are determined 
through commutation relations with the Liouvillian. We 
have found the steady state and the expansion coefficients 
by using step-down superoperators and, the remaining 
normal modes, by stepping the steady state up. 
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APPENDIX A: COMPLETENESS OF THE BASIS 
FORMED BY EIGENSTATES OF ONE 
NON-THERMAL MODE 

We want to demonstrate that the eigenstates i?™>" 
with positive integer indices form a complete set, i.e., 
they can expand any target operator p. For this demon- 
stration we will use a straightforward calculation, i.e., we 
will use the explicit form of these eigenstates and of the 
expansion coefficients and suitably sum the series. 
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The expansion coefficients can be developed by per- 
forming the trace of the target operator p in the basis 
of the oscillator eigenstates and applying the operators a 
and presents in the step-down superoperators (11) 



tr| — ^= p 



'n\ vm\ 



1=0 



{m + l)\{n + l)\ 



ml n! 
x{m + l\p\n + l) 



Using the explicit form of the eigenstates R" 
Eq. (23), we can sum the series 



(Al) 



given m 



oo oo 

EE 

m—0 n— 



tr ■ 



p]R" 



OO OO OO OO 

-EEEE 



(-1)'= /(m + /)! {n + iy. 



kill 



k—O m—k n—k / — 

X (to + l\p\n + 1} \m — k){n — k\ , 



TO — k)l {n — k)l 



(A2) 



where we have suitably reordered the summations. By 
changing the sum indices m^m + k, n^n + k and 
I ^ I — k and reordering the summations, we get 



EE*^ 

m=0 n=0 

oo oo oo ? 



Ml 



n\ VTO! 



pi?" 



kl{l-k)\ 



m=0 n=0 1=0 k=0 

X {m + l\p\n + l) \m){n 



{-If {m + iy.{n + l)l 



ml 



(A3) 



The summation over k is easily done resulting in ^;,0) 
which eliminates the summation over L It is easy to rec- 
ognize that the resulting expression is exactly the ex- 
pansion of the operator p into the basis of the oscillator 
eigenstates. This allows us to obtain 



CO oo 

EE 

m=0 n=0 



tr 



A/"!! M[ 



nl VTO! 



(A4) 



This last expression shows that the set of the eigenstates 
^m,n positive integer indices is complete. 

In this proof, it was essential to reorder the summa- 
tions, but this is not always possible because it requires 
certain conditions. Then the above demonstration will 
not be always valid. As criterion to determine the valid- 
ity of (A4), we have that this equality will occur if and 
only if each matrix element of the expansion defined in its 
left-hand side converges. A simple and concrete example 
of this fact appears in the case of a thermal distribution 



°° °° / A/"" A^' 

E E^M ^^=f /'(^^ 

m=0 n=0 



oo oo 



EE(-i) 

fe=0 n=k 



n—k 



n\ vm! 
nl 



kl{n-k)l 



n''\k){k\. 



(A6) 



This expansion correctly converges for (A5) when n < 1, 
due to identity (48), but diverges when n > 1. 

The completeness for the remaining cases can be also 
demonstrated, however it is more arduous and therefore 
will be omitted. 
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= E 



-\k){k\. 



Its expansion into the eigenstates will give us 



(A5) 
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